## Chapter 6 (Evaluations) 
## Spring 2022 

################
## Figure 6.1 ##
################

jeps <- read_dta(file="~/Dropbox/Affective Polarization/COVID/JEPS/Replication Data/Aff. Pol. Exp. Data.dta") 

## Code up partisan animosity scale & experimental DVs 
jeps <- jeps %>% 
  mutate(allout2 = (feelingtherm_outgroup +patriotic_outgroup +intelligent_outgroup +honest_outgroup +
                      openminded_outgroup +generous_outgroup +hypocritical_outgroup + selfish_outgroup +
                      mean_outgroup +trust_outgroup +friends_outgroup +neighbors_outgroup +marry_outgroup)/13, 
         negparty = 1 - allout2, 
         rev_prep = 5 - expprep, 
         rev_future = 5 - expfuture) %>% 
  rowwise() %>% 
  mutate(exp_dv = mean(c(expconf,rev_prep,rev_future), na.rm=T)) %>% 
  ungroup() 

## quantiles for x-axis 
jeps_quantiles <- quantile(jeps$negparty, probs=c(0.05,0.10,0.33,0.5,0.67,0.90,0.95), na.rm=T)  

## create a variable that's party x condition 
jeps <- jeps %>% 
  mutate(trump_exp = case_when( 
    exptrump == 1 & dem == 1 ~ 1, 
    exptrump == 0 & dem == 1 ~ 2, 
    exptrump == 1 & dem == 0 ~ 3, 
    exptrump == 0 & dem == 0 ~ 4))   

## create factor for US vs. Trump 
jeps$trump_exp_factor <- factor(jeps$exptrump,
                                labels=c("United States","Trump"))

pdf(file = "Figures/Druckman_etal_Figure_6_1.pdf", 
    height = 8.5,
    width = 11) 
print(jeps %>% 
        ## remove pure Independents 
        filter(pid != 4) %>% 
        ## filter out the top/bottom 5% to match other graphs 
        filter(negparty > ap_quantiles[1] & negparty < ap_quantiles[7]) %>%    
        ## remove those who didn't see either version (not in this wave) 
        filter(!is.na(trump_exp_factor)) %>% 
        ## draw the graph! 
        ggplot(.) + 
        geom_smooth(se = T, method = "gam",
                    aes(x = negparty, y = exp_dv, 
                        color = factor(rep),
                        fill = factor(rep)))+  
        scale_color_manual(name = "", 
                           labels = c("Democrat","Republican"), 
                           values = c("black","grey70")) + 
        scale_fill_manual(name = "", 
                           labels = c("Democrat","Republican"), 
                           values = c("black","grey70")) + 
        scale_x_continuous(breaks = jeps_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        xlab("Animus") +
        ylab("Approval of Pandemic Response") + 
        ylim(c(1,4)) + 
        #guides(linetype = guide_legend(override.aes = list(color="black"))) +  
        guides(color = guide_legend(),
               fill = guide_legend()) + 
        facet_grid(cols=vars(trump_exp_factor)) + 
        #ggtitle("Partisan Animosity and Approval of Handling COVID-19") + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none",
              panel.grid = element_blank())) 
dev.off() 

# ## Here's the graph with the loess instead of OLS for the appendix: 
# pdf(file="Figures/Extra/jeps_with_loess.pdf",
#     height = 11,
#     width = 8.5)
# print(jeps %>% 
# ## remove pure Independents 
# filter(pid != 4) %>% 
#   filter(negparty > 0.265) %>% 
# ## draw the graph! 
# ggplot(.) + 
# geom_smooth(se = F, 
# aes(x = negparty, y = exp_dv, 
# color = factor(dem, 
# labels = c("Republicans","Democrats")), 
# linetype = factor(exptrump, 
# labels = c("U.S.","Trump")))) + 
# scale_color_manual(name = "", 
# labels = c("Republican","Democrat"), 
# values = c("grey60","black")) + 
# scale_linetype_manual(name = "",
# labels = c("U.S.","Trump"),
# values = c("dashed","solid")) +  
# scale_x_continuous(breaks = jeps_quantiles, 
# labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
# xlab("Partisan Animus") +
# ylab("Approval of Pandemic Response") + 
# ylim(c(1,4)) + 
# guides(linetype = guide_legend(override.aes = list(color="black"))) +  
# ggtitle("Partisan Animosity & Approval of Trump vs. U.S.") + 
# theme_bw() + 
# theme(plot.title = element_text(hjust=0.5),
# legend.position = "bottom",
# panel.grid = element_blank())) 
# dev.off() 

############################################
## Figure 6.2 & 6.3: Biden/Trump Approval ## 
############################################

## create an indicator for party*wave 
data <- data %>% 
  mutate(party_year = case_when(
    dem == 1 & year == 2 ~ 1, 
    dem == 1 & year == 3 ~ 2, 
    dem == 1 & year == 4 ~ 3, 
    dem == 0 & year == 2 ~ 4, 
    dem == 0 & year == 3 ~ 5, 
    dem == 0 & year == 4 ~ 6))


pdf(file="Figures/Druckman_etal_figure_6_2.pdf", 
    height = 8.5, 
    width = 11) 
## Trump Approval, Faceted 
print(data %>% 
  #filter(partyid_w1 != 4 & partychange == 0) %>% 
  filter(year > 1) %>% 
  ggplot(aes(apstable, trumpapr,
             group = rep)) + 
  geom_smooth(aes(color = as.factor(rep),
                  fill = as.factor(rep)), se = T) + 
  scale_color_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) + 
  scale_fill_manual(name = "",
                       values = c("black","grey70"), 
                       labels = c("Democrat","Republican")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                "Neither Approve \n nor Disapprove",
                                "Somewhat \n Approve", "Strongly \n Approve")) + 
  xlab("Animus") +
  ylab("Trump Approval") + 
  #ylim(c(1,5)) + 
  #ggtitle("Approval of Trump's Handling of COVID-19") + 
  labs(color = "") + 
  theme_bw() + 
  facet_grid(cols=vars(factor_year))+
  guides(color = guide_legend(nrow=2,byrow=T),
         linetype = guide_legend(nrow=2,byrow=T)) + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank()))  
dev.off()

## Biden approval, faceted 
pdf(file="Figures/Druckman_etal_figure_6_3.pdf", 
    height = 8.5, 
    width = 11) 
## Biden approval on asked in April 2020/April 2021 
print(data %>% 
        #filter(partyid_w1 != 4 & partychange == 0) %>% 
        filter(year != 1 & year !=3) %>% 
        ggplot(aes(apstable, bidenapr,
                   group = rep)) + 
        geom_smooth(aes(color = as.factor(rep),
                        fill = as.factor(rep)), se = T) + 
        scale_color_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Democrat","Republican")) + 
        scale_fill_manual(name = "",
                          values = c("black","grey70"), 
                          labels = c("Democrat","Republican")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        scale_y_continuous(breaks = c(1:5),
                           limits = c(1,5),
                           labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                      "Neither Approve \n nor Disapprove",
                                      "Somewhat \n Approve", "Strongly \n Approve")) + 
        xlab("Animus") +
        ylab("Biden Approval") + 
        #ylim(c(1,5)) + 
        #ggtitle("Approval of Biden's Handling of COVID-19") + 
        labs(color = "") + 
        theme_bw() + 
        facet_grid(cols=vars(factor_year))+
        guides(color = guide_legend(nrow=2,byrow=T),
               linetype = guide_legend(nrow=2,byrow=T)) + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid = element_blank()))  
dev.off()

#######################################################
## Figure 6.4 & 6.5: State & Local Officials, Mayors ## 
#######################################################

## State and Local Officials 
pdf(file="Figures/Druckman_etal_Figure_6_4.pdf", 
    height = 11, 
    width = 8.5) 
print(data %>% 
  filter(year == 2) %>% 
  ggplot(aes(apstable, cvastate,
             group = dem)) + 
  geom_smooth(aes(color = as.factor(dem),
                  fill = as.factor(dem)), se = T) + 
  scale_color_manual(values = c("grey70","black"), 
                     labels = c("Republicans","Democrats")) + 
  scale_fill_manual(values = c("grey70","black"), 
                     labels = c("Republicans","Democrats")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                "Neither Approve \n nor Disapprove",
                                "Somewhat \n Approve", "Strongly \n Approve")) + 
  xlab("Animus") +
  ylab("State and Local Officials Approval") + 
  #ylim(c(1,5)) + 
  #ggtitle("Approval of State and Local Officials Handling of COVID-19 (April 2020)") + 
  labs(color = "") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank())) 
dev.off() 

## Mayoral Approval, Faceted 
pdf(file="Figures/Druckman_etal_figure_6_5.pdf", 
    height = 8.5, 
    width = 11) 
print(data %>% 
  filter(year > 2) %>% 
  ggplot(aes(apstable, mayorapr,
             group = rep)) + 
  geom_smooth(aes(color = as.factor(rep),
                  fill = as.factor(rep)), se = T) + 
  scale_color_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) + 
  scale_fill_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                "Neither Approve \n nor Disapprove",
                                "Somewhat \n Approve", "Strongly \n Approve")) + 
  xlab("Animus") +
  ylab("Mayoral Approval") + 
  #ylim(c(1,5)) + 
  #ggtitle("Approval of Your Mayor's Handling of COVID-19") + 
  labs(color = "") + 
  facet_grid(cols=vars(factor_year)) + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank()))  
dev.off()

################
## Figure 6.6 ## 
################

pdf(file = "Figures/Druckman_etal_Figure_6_6.pdf",
    height = 11,
    width = 8.5) 
print(data %>% 
  filter(year == 2) %>% 
  ggplot(aes(apstable, cvacongress,
             group = dem)) + 
  geom_smooth(aes(color = as.factor(dem),
                  fill = as.factor(dem)), se = T) + 
  scale_color_manual(values = c("grey70","black"), 
                     labels = c("Republicans","Democrats")) + 
  scale_fill_manual(values = c("grey70","black"), 
                       labels = c("Republicans","Democrats")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                "Neither Approve \n nor Disapprove",
                                "Somewhat \n Approve", "Strongly \n Approve")) + 
  xlab("Animus") +
  ylab("Members of Congress Approval") + 
  #ylim(c(1,5)) + 
  #ggtitle("Approval of Congress's Handling of COVID-19") + 
  labs(color = "") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank()))  
dev.off() 

################################
## Figure 6.7: Vaccine Credit ## 
################################

## combine distribution and development of the vaccine (scale strongly) 
psych::alpha(cbind(data$w7outpresdevvaccine,data$w7outpresdistrvaccine))
psych::alpha(cbind(data$w7inpresdevvaccine,data$w7inpresdistrvaccine))
psych::alpha(cbind(data$w7govdevvaccine,data$w7govdistrvaccine))

## for govenor: read in excel file w/ partisanship (to make same-party/other-party labels)
gov <- read_xlsx(path = "~/Dropbox/Affective Polarization/Book Analyses/OLD MATERIALS/Chapter 6 Evaluations/gov_partisanship.xlsx") %>% clean_names(.)
data_with_govparty <- left_join(x = data, 
                                y = gov, 
                                by = c("statefips" = "state_fips")) 
## any fips codes botched? 
anti_join(x = gov, 
          y = data, 
          by = c("state_fips" = "statefips")) 
## all good 


## Make as a Faceted Graph 
## First step: reshape presidential/governor credit/blame for vaccine wide --> long 
## Second step: do the graph! 
pdf(file="Figures/Druckman_etal_figure_6_7.pdf", 
    height = 11, 
    width = 8.5) 
print(data_with_govparty %>% 
  filter(partyid != 4) %>% 
  filter(year == 4) %>% 
  rowwise(.) %>% 
  mutate(vax_gov = mean(c(w7govdevvaccine,w7govdistrvaccine),na.rm=T),
         vax_pres_inpty = mean(c(w7inpresdevvaccine,w7inpresdistrvaccine),na.rm=T), 
         vax_pres_outpty = mean(c(w7outpresdevvaccine,w7outpresdistrvaccine),na.rm=T)) %>%  
  ungroup(.) %>% 
  mutate(party_match = ifelse((gov_party == "D" & dem == 1) |
                                (gov_party == "R" & dem == 0),1,0),
         vax_gov_inpty = ifelse(party_match == 1, vax_gov, NA),
         vax_gov_outpty = ifelse(party_match == 0, vax_gov, NA)) %>% 
  select(c(apstable,vax_pres_inpty,vax_pres_outpty,vax_gov_inpty,vax_gov_outpty)) %>% 
  pivot_longer(cols = c(vax_pres_inpty,vax_pres_outpty,vax_gov_inpty,vax_gov_outpty),
               names_to = "group",
               values_to = "eval") %>% 
  separate(group, c("vax","leader","party"), sep="_") %>% 
  select(-vax) %>% 
  mutate(leader_factor = factor(leader,
                                levels = c("pres","gov"),
                                labels = c("President","Governor"))) %>% 
  ## now graph 
  ggplot(aes(x = apstable, y = eval, group = party)) + 
  geom_smooth(aes(linetype = factor(party)),
              col="black",se=T) + 
  scale_linetype_manual(values = c("solid","dashed"), 
                        labels = c("In-Party","Out-Party")) +
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Definitely \n Hurt","Probably \n Hurt",
                                "Neither Helped \n nor Hurt",
                                "Probably \n Helped", "Definitely \n Helped")) + 
  xlab("Animus") +
  ylab("President/Governor Helped or Hurt COVID-19 Vaccines") + 
  #ylim(c(1,5)) + 
  #ggtitle("Help or Hurt the Vaccine Program") + 
  facet_grid(cols = vars(leader_factor)) + 
  labs(lty = "") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank())) 
dev.off()
  
################################
## Figure 6.8: Fauci Approval ## 
################################

fauci <- data %>% 
        filter(year > 1) %>% 
        #filter(partyid_w1 != 4 & partychange == 0) %>% 
        ggplot(aes(apstable, fauciapr,
                   group = rep)) + 
        geom_smooth(aes(color = as.factor(rep),
                        fill = as.factor(rep)), se = T) + 
        scale_color_manual(name = "",
                           values = c("black","grey70"), 
                           labels = c("Democrat","Republican")) + 
        scale_fill_manual(name = "",
                          values = c("black","grey70"), 
                          labels = c("Democrat","Republican")) + 
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        scale_y_continuous(breaks = c(1:5),
                           limits = c(1,5),
                           labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                      "Neither Approve \n nor Disapprove",
                                      "Somewhat \n Approve", "Strongly \n Approve")) + 
        xlab("Animus") +
        ylab("Fauci Approval") + 
        #ylim(c(1,5)) + 
        #ggtitle("Approval of Dr. Fauci's Handling of COVID-19") + 
        labs(color = "") + 
        #guides(color = guide_legend(nrow=2,byrow=T),
        #       linetype = guide_legend(nrow=2,byrow=T)) + 
        facet_grid(cols = vars(factor_year)) + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid = element_blank()) 

## Graph 
pdf(file="Figures/Druckman_etal_Figure_6_8.pdf", 
    height = 8.5,
    width = 11) 
print(fauci)
dev.off() 

## look at data to understand gaps over time 
## these give you the max gap, by panel 
fauci_data <- ggplot_build(fauci)$data[[1]] 
fauci_data$y[80] - fauci_data$y[160]  ## 0.34 (April 2020, this is panel 1) 
fauci_data$y[240] - fauci_data$y[320]  ## 1.7 (October 2020, this is panel 2) 
fauci_data$y[400] - fauci_data$y[480]  ## 2.5 (April 2021, this is panel 3) 

## gap is 0.34 in April 2020, where is that in October? 
quantile(data$apstable, probs=(0.15))  
fauci_data$y[175] - fauci_data$y[255]  ## 0.38, about as close as you can get 

#########################################
## Figure 6.9 & 6.10: CDC/FDA Approval ##
#########################################

## CDC as a Facet
pdf(file="Figures/Druckman_etal_Figure_6_9.pdf", 
    height = 8.5, 
    width = 11) 
print(data %>% 
  filter(year > 1) %>% 
  #filter(partyid_w1 != 4 & partychange == 0) %>% 
  ggplot(aes(apstable, cdcapr,
             group = rep)) + 
  geom_smooth(aes(color = as.factor(rep),
                  fill = as.factor(rep)), se = T) + 
  scale_color_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) + 
  scale_fill_manual(name = "",
                    values = c("black","grey70"), 
                    labels = c("Democrat","Republican")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                "Neither Approve \n nor Disapprove",
                                "Somewhat \n Approve", "Strongly \n Approve")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  xlab("Animus") +
  ylab("CDC Approval") + 
  #ylim(c(1,5)) + 
  #ggtitle("Approval of The CDC's Handling of COVID-19") + 
  labs(color = "") + 
  facet_grid(cols = vars(factor_year)) + 
  theme_bw() + 
  guides(color = guide_legend(nrow=2,byrow=T),
         linetype = guide_legend(nrow=2,byrow=T)) + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank())) 
dev.off() 

## FDA as a Facet 
## FDA Approval (only asked in Waves 3-4)
pdf(file="Figures/Druckman_etal_figure_6_10.pdf", 
    height = 8.5, 
    width = 11) 
print(data %>% 
  filter(year > 2) %>% 
  ggplot(aes(apstable, fdaapr,
             group = rep)) + 
  geom_smooth(aes(color = as.factor(rep),
                  fill = as.factor(rep)), se = T) + 
  scale_color_manual(name = "",
                     values = c("black","grey70"), 
                     labels = c("Democrat","Republican")) + 
  scale_fill_manual(name = "",
                    values = c("black","grey70"), 
                    labels = c("Democrat","Republican")) + 
  scale_x_continuous(breaks = ap_quantiles, 
                     labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
  scale_y_continuous(breaks = c(1:5),
                     limits = c(1,5),
                     labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                "Neither Approve \n nor Disapprove",
                                "Somewhat \n Approve", "Strongly \n Approve")) + 
  xlab("Animus") +
  ylab("FDA Approval") + 
  #ylim(c(1,5)) + 
  #ggtitle("Approval of The FDA's Handling of COVID-19") + 
  labs(color = "") + 
  theme_bw() + 
  facet_grid(cols = vars(factor_year)) + 
  theme(plot.title = element_text(hjust=0.5),
        legend.position = "none", 
        panel.grid = element_blank()))  
dev.off() 


###########################################################
## Figure 6.11: Severity and Scientific Advisor Approval ## 
###########################################################

## group severity into two bins, too much jumpiness with all four categories 
data$cm <- ifelse(data$covidmax == 1 | data$covidmax == 2,1,
                  ifelse(data$covidmax==3 | data$covidmax == 4,2,NA)) 
## put everyone into one plot 
data <- data %>% 
  mutate(cm_party = case_when( 
    cm == 1 & dem == 0 ~ 1, 
    cm == 2 & dem == 0 ~ 2, 
    cm == 1 & dem == 1 ~ 3, 
    cm == 2 & dem == 1 ~ 4)) 

## Pool October/April 21 waves to make it clearer (less bumpy)
## first make a factor variable for the labels 
data$sev_label <- factor(data$cm, 
                         labels = c("Low Consequences",
                                    "High Consequences"))


## do all three on one plot 
pdf(file="Figures/Druckman_etal_Figure_6_11.pdf",
    height = 8.5,
    width = 11)
print(data %>% 
        filter(year > 2 & !is.na(sev_label)) %>% 
        select(c(apstable,fauciapr,cdcapr,fdaapr,dem,sev_label)) %>% 
        pivot_longer(cols = c(fauciapr,cdcapr,fdaapr),
                     names_to = "group",
                     values_to = "eval") %>% 
        mutate(leader_factor = factor(group, 
                                      levels = c("fauciapr","cdcapr","fdaapr"),
                                      labels = c("Dr. Fauci","CDC","FDA"))) %>% 
        ## now graph 
        ggplot(aes(x = apstable, y = eval, group = dem)) + 
        geom_smooth(aes(color = as.factor(dem),
                        fill = as.factor(dem)),se=T) + 
        scale_color_manual(values = c("grey70","black"), 
                              labels = c("Republican","Democrat")) +
        scale_fill_manual(values = c("grey70","black"), 
                           labels = c("Republican","Democrat")) +
        scale_x_continuous(breaks = ap_quantiles, 
                           labels = c("5th \n pct","10th \n pct","33rd \n pct","50th \n pct","67th \n pct", "90th \n pct","95th \n pct")) + 
        scale_y_continuous(breaks = c(1:5),
                            limits = c(1,5),
                            labels = c("Strongly \n Disapprove","Somewhat \n Disapprove",
                                       "Neither Approve \n nor Disapprove",
                                       "Somewhat \n Approve", "Strongly \n Approve")) + 
        xlab("Animus") +
        ylab("Approval") + 
        #ylim(c(1,5)) + 
        #ggtitle("Salience and Scientific Advisor Approval") + 
        facet_grid(cols = vars(sev_label),
                   rows = vars(leader_factor)) + 
        labs(lty = "") + 
        theme_bw() + 
        theme(plot.title = element_text(hjust=0.5),
              legend.position = "none", 
              panel.grid = element_blank())) 
dev.off()

########################################
## Panel Regressions for the Appendix ##
######################################## 

m1 <- felm(trumpapr ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
m2 <- felm(bidenapr ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
m3 <- felm(mayorapr ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
m4 <- felm(fauciapr ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
m5 <- felm(cdcapr ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)
m6 <- felm(fdaapr ~ year*apstable*rep | uniqueid + year, 
           subset = partyid_w1 != 4 & partychange == 0,
           data = data)

## output to tables 
## output this to a table 
models <- list()
models[["Trump"]] <- m1  
models[["Biden"]] <- m2 
models[["Mayor"]] <- m3
models[["Fauci"]] <- m4
models[["CDC"]] <- m5
models[["FDA"]] <- m6

## add the N 
rows <- tribble(~term, ~"m1",~"m2",~"m3",~"m4",~"m5",~"m6",
                'N',m1$N,m2$N,m3$N,m4$N,m5$N,m6$N)
attr(rows,'position') <- 15 

## clean up the GOF output (just N & R-Squared)
gm <- modelsummary::gof_map 
gm$omit <- TRUE 
gm$omit[5] <- FALSE 

msummary(models,
         add_rows = rows, 
         coef_map = c("apstable" = "Animus",
                      "rep" = "Republican",
                      "year" = "Wave",
                      "year:rep" = "Animosity*Republican",
                      "year:apstable" = "Animosity*Wave",
                      "apstable:rep" = "Republican*Wave", 
                      "year:apstable:rep" = "Animosity*Republican*Wave"), 
         gof_map = gm, 
         stars = T, 
         output =  "Tables/raw/table_A6_1.html")   


## Idea: look at changes in COVID severity, see if it relates to changes in Fauci approval

## First, calculate the change in the DVs and the covidmax variable 
data <- data %>% 
  group_by(respondent_id) %>% 
  mutate(cm_lag = dplyr::lag(covidmax, n = 1, order_by = year), 
         fauci_lag = dplyr::lag(fauciapr, n = 1, order_by = year), 
         cdc_lag = dplyr::lag(cdcapr, n = 1, order_by = year),
         fda_lag = dplyr::lag(fdaapr, n = 1 , order_by = year)) %>% 
  rowwise(.) %>% 
  mutate(delta_cm = covidmax - cm_lag, 
         delta_fauci = fauciapr - fauci_lag, 
         delta_cdc = cdcapr - cdc_lag, 
         delta_fda = fdaapr - fda_lag) %>% 
  ungroup(.) 

m1 <- lm(delta_fauci ~ delta_cm, data = data) 
m2 <- lm(delta_cdc ~ delta_cm, data = data)
m3 <- lm(delta_fda ~ delta_cm, data = data) 
m4 <- lm(delta_fauci ~ delta_cm, subset = rep == 1, data = data) 
m5 <- lm(delta_cdc ~ delta_cm, subset = rep == 1, data = data)
m6 <- lm(delta_fda ~ delta_cm, subset = rep == 1, data = data) 
## overall effects for everyone, but among Rs only CDC effects (Fauci p = 0.15)

## export to a table 
models <- list()
models[["Dr. Fauci \n (All)"]] <- m1 
models[["CDC \n (All)"]] <- m2 
models[["FDA \n (All)"]] <- m3 
models[["Fauci \n (Rs Only)"]] <- m4 
models[["CDC \n (Rs Only)"]] <- m5
models[["FDA \n (Rs Only)"]] <- m6 

## clean up the GOF statistics (just add N/R2)
gm <- gof_map
gm$omit <- T 
gm$omit[2] <- F
gm$omit[5] <- F 

msummary(models,
         coef_map = c("(Intercept)" = "Intercept",
                      "delta_cm" = "Change in COVID Severity"), 
         stars = T, 
         #gof_map = gm, 
         #escape = F, 
         output =  "Tables/raw/Table_A6_2.html")
